
#define LIMIT -999


#include <stdlib.h>

#include <stdio.h>
#include <string.h>
#include <math.h>


#include <sys/time.h>
#include "mytime.h"


//#include <omp.h>


#define OPENMP

#define NUM_THREAD 4

////////////////////////////////////////////////////////////////////////////////
// declaration, forward
//void runTest( int argc, char** argv);
void runTest_llc( int max_rows, int max_cols, int penalty, int * input_itemsets, int *referrence);
void runTest_omp( int max_rows, int max_cols, int penalty, int * input_itemsets, int *referrence);
void prepare (int * input_itemsets, int * b, int max_rows, int max_cols);
void usage (int argc, char**argv);


int check_matrix (int *a, int *b, int size) {
   int i = 0;
   for (i; i < size; i++) 
	if (a[i] != b[i]) {
		printf ("Error!!! Llc matrix is not equal to OMP matrix. LLC [%d] - %d \t OMP [%d] - %d\n", i, a[i], i, b[i]);
		return -1;
	}
   return 0;
}


int maximum( int a,
		 int b,
		 int c){

	int k;
	if( a <= b )
		k = b;
	else 
	k = a;

	if( k <=c )
	return(c);
	else
	return(k);
}


int blosum62[24][24] = {
{ 4, -1, -2, -2,  0, -1, -1,  0, -2, -1, -1, -1, -1, -2, -1,  1,  0, -3, -2,  0, -2, -1,  0, -4},
{-1,  5,  0, -2, -3,  1,  0, -2,  0, -3, -2,  2, -1, -3, -2, -1, -1, -3, -2, -3, -1,  0, -1, -4},
{-2,  0,  6,  1, -3,  0,  0,  0,  1, -3, -3,  0, -2, -3, -2,  1,  0, -4, -2, -3,  3,  0, -1, -4},
{-2, -2,  1,  6, -3,  0,  2, -1, -1, -3, -4, -1, -3, -3, -1,  0, -1, -4, -3, -3,  4,  1, -1, -4},
{ 0, -3, -3, -3,  9, -3, -4, -3, -3, -1, -1, -3, -1, -2, -3, -1, -1, -2, -2, -1, -3, -3, -2, -4},
{-1,  1,  0,  0, -3,  5,  2, -2,  0, -3, -2,  1,  0, -3, -1,  0, -1, -2, -1, -2,  0,  3, -1, -4},
{-1,  0,  0,  2, -4,  2,  5, -2,  0, -3, -3,  1, -2, -3, -1,  0, -1, -3, -2, -2,  1,  4, -1, -4},
{ 0, -2,  0, -1, -3, -2, -2,  6, -2, -4, -4, -2, -3, -3, -2,  0, -2, -2, -3, -3, -1, -2, -1, -4},
{-2,  0,  1, -1, -3,  0,  0, -2,  8, -3, -3, -1, -2, -1, -2, -1, -2, -2,  2, -3,  0,  0, -1, -4},
{-1, -3, -3, -3, -1, -3, -3, -4, -3,  4,  2, -3,  1,  0, -3, -2, -1, -3, -1,  3, -3, -3, -1, -4},
{-1, -2, -3, -4, -1, -2, -3, -4, -3,  2,  4, -2,  2,  0, -3, -2, -1, -2, -1,  1, -4, -3, -1, -4},
{-1,  2,  0, -1, -3,  1,  1, -2, -1, -3, -2,  5, -1, -3, -1,  0, -1, -3, -2, -2,  0,  1, -1, -4},
{-1, -1, -2, -3, -1,  0, -2, -3, -2,  1,  2, -1,  5,  0, -2, -1, -1, -1, -1,  1, -3, -1, -1, -4},
{-2, -3, -3, -3, -2, -3, -3, -3, -1,  0,  0, -3,  0,  6, -4, -2, -2,  1,  3, -1, -3, -3, -1, -4},
{-1, -2, -2, -1, -3, -1, -1, -2, -2, -3, -3, -1, -2, -4,  7, -1, -1, -4, -3, -2, -2, -1, -2, -4},
{ 1, -1,  1,  0, -1,  0,  0,  0, -1, -2, -2,  0, -1, -2, -1,  4,  1, -3, -2, -2,  0,  0,  0, -4},
{ 0, -1,  0, -1, -1, -1, -1, -2, -2, -1, -1, -1, -1, -2, -1,  1,  5, -2, -2,  0, -1, -1,  0, -4},
{-3, -3, -4, -4, -2, -2, -3, -2, -2, -3, -2, -3, -1,  1, -4, -3, -2, 11,  2, -3, -4, -3, -2, -4},
{-2, -2, -2, -3, -2, -1, -2, -3,  2, -1, -1, -2, -1,  3, -3, -2, -2,  2,  7, -1, -3, -2, -1, -4},
{ 0, -3, -3, -3, -1, -2, -2, -3, -3,  3,  1, -2,  1, -1, -2, -2,  0, -3, -1,  4, -3, -2, -1, -4},
{-2, -1,  3,  4, -3,  0,  1, -1,  0, -3, -4,  0, -3, -3, -2,  0, -1, -4, -3, -3,  4,  1, -1, -4},
{-1,  0,  0,  1, -3,  3,  4, -2,  0, -3, -3,  1, -1, -3, -1,  0, -1, -3, -2, -2,  1,  4, -1, -4},
{ 0, -1, -1, -1, -2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -2,  0,  0, -2, -1, -1, -1, -1, -1, -4},
{-4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4, -4,  1}
};

double gettime() {
  struct timeval t;
  gettimeofday(&t,NULL);
  return t.tv_sec+t.tv_usec*1e-6;
}

////////////////////////////////////////////////////////////////////////////////
// Program main
////////////////////////////////////////////////////////////////////////////////
double ex_time;
double omp_time;

int max_rows;
int max_cols;
int penalty;
int * matrix_llc;
int * matrix_omp;
int * referrence;



int
main( int argc, char** argv) 
{

    // the lengths of the two sequences should be able to divided by 16.
    // And at current stage  max_rows needs to equal max_cols
    if (argc == 3)
    {
	max_rows = atoi(argv[1]);
	max_cols = atoi(argv[1]);
	penalty = atoi(argv[2]);
	//omp_num_threads = atoi(argv[3]);
    }
    else{
	usage(argc, argv);
    }

    max_rows = max_rows + 1;
    max_cols = max_cols + 1;
    referrence = (int *)malloc( max_rows * max_cols * sizeof(int) );
    matrix_llc = (int *)malloc( max_rows * max_cols * sizeof(int) );
    matrix_omp = (int *)malloc( max_rows * max_cols * sizeof(int) );
	

    if (!matrix_omp)
	fprintf(stderr, "error: can not allocate memory");
    prepare(matrix_llc, matrix_omp, max_rows, max_cols);
    if (check_matrix(matrix_llc, matrix_omp, max_rows*max_cols) == 0)
	printf ("The initial matrix are the same\n");
    runTest_llc(max_rows, max_cols, penalty, matrix_llc, referrence); 
    printf("LLC Execution time: %g (s)\n", ex_time);
    runTest_omp(max_rows, max_cols, penalty, matrix_omp, referrence); 
    printf("OMP Execution time: %g (s)\n", omp_time);
    if (check_matrix(matrix_llc, matrix_omp, max_rows*max_cols) == 0)
	printf ("The result have been verified\n");

    return EXIT_SUCCESS;
}

void usage(int argc, char **argv)
{
	fprintf(stderr, "Usage: %s <max_rows/max_cols> <penalty> <num_threads>\n", argv[0]);
	fprintf(stderr, "\t<dimension>      - x and y dimensions\n");
	fprintf(stderr, "\t<penalty>        - penalty(positive integer)\n");
	//fprintf(stderr, "\t<num_threads>    - no. of threads\n");
	exit(1);
}


void prepare (int * input_itemsets, int * b, int max_rows, int max_cols)
{


    int i,j;
    srand ( 7 );

    for ( i = 0 ; i < max_cols; i++){
		for ( j = 0 ; j < max_rows; j++){
			input_itemsets[i*max_cols+j] = 0;
		}
	}

	printf("Start Needleman-Wunsch\n");

	for(  i=1; i< max_rows ; i++){    //please define your own sequence. 
       input_itemsets[i*max_cols] = rand() % 10 + 1;
	}
    for(  j=1; j< max_cols ; j++){    //please define your own sequence.
       input_itemsets[j] = rand() % 10 + 1;
	}


	for ( i = 1 ; i < max_cols; i++){
		for ( j = 1 ; j < max_rows; j++){
		referrence[i*max_cols+j] = blosum62[input_itemsets[i*max_cols]][input_itemsets[j]];
		}
	}

    for(  i = 1; i< max_rows ; i++)
       input_itemsets[i*max_cols] = -i * penalty;
	for( j = 1; j< max_cols ; j++)
       input_itemsets[j] = -j * penalty;
    for (i = 0; i < max_rows*max_cols;i++)
        b[i] = input_itemsets[i];
}


////////////////////////////////////////////////////////////////////////////////
//! Run a simple test for CUDA
////////////////////////////////////////////////////////////////////////////////
void
runTest_llc( int max_rows, int max_cols, int penalty, int * input_itemsets, int *referrence) 
{
    int idx, index;
    //int max_rows, max_cols, penalty,idx, index;
    //int *input_itemsets, *output_itemsets, *referrence;
    //int *matrix_cuda, *matrix_cuda_out, *referrence_cuda;
    int size;
    int omp_num_threads;
    int i, j;
	
    



	
	//Compute top-left matrix 
	//printf("Num of threads: %d\n", omp_num_threads);
	printf("Processing top-left matrix\n");
    CLOCK_TYPE chrono;
    CLOCK_Start(chrono);

    #pragma llc context name("compute_main") copy_in(input_itemsets[max_rows * max_cols], penalty, i, index, max_cols, referrence[max_rows * max_cols], idx) copy_out(input_itemsets[max_rows * max_cols], referrence[max_rows * max_cols]) copy_out(input_itemsets[max_rows * max_cols])
    {
    for(i = 0 ; i < max_cols-2 ; i++){
       #pragma llc for shared(input_itemsets, referrence, i,max_cols,penalty) private(idx, index) 
		for( idx = 0 ; idx <= i ; idx++){
		 index = (idx + 1) * max_cols + (i + 1 - idx);
         input_itemsets[index]= maximum( input_itemsets[index-1-max_cols]+ referrence[index], 
			                             input_itemsets[index-1]         - penalty, 
									     input_itemsets[index-max_cols]  - penalty);

		}
	}
	printf("Processing bottom-right matrix\n");
    //Compute bottom-right matrix 
	for(i = max_cols - 4 ; i >= 0 ; i--){
       #pragma llc for shared(input_itemsets, referrence, i,max_cols,penalty) private(idx, index) 
       for( idx = 0 ; idx <= i ; idx++){
	      index =  ( max_cols - idx - 2 ) * max_cols + idx + max_cols - i - 2 ;
		  input_itemsets[index]= maximum( input_itemsets[index-1-max_cols]+ referrence[index], 
			                              input_itemsets[index-1]         - penalty, 
									      input_itemsets[index-max_cols]  - penalty);
	      }

	}
    }
    CLOCK_End(chrono, ex_time);
   	

}


////////////////////////////////////////////////////////////////////////////////
//! Run a simple test for OMP
////////////////////////////////////////////////////////////////////////////////
void
runTest_omp( int max_rows, int max_cols, int penalty, int * input_itemsets, int *referrence) 
{
    int idx, index;
    //int max_rows, max_cols, penalty,idx, index;
    //int *input_itemsets, *output_itemsets, *referrence;
    //int *matrix_cuda, *matrix_cuda_out, *referrence_cuda;
    int size;
    int omp_num_threads;
    int i, j;
	
    



	
	//Compute top-left matrix 
	//printf("Num of threads: %d\n", omp_num_threads);
	printf("Processing top-left matrix\n");
    CLOCK_TYPE chrono;
    CLOCK_Start(chrono);

    {
    for(i = 0 ; i < max_cols-2 ; i++){
       #pragma omp parallel for shared(input_itemsets, referrence, i,max_cols,penalty) private(idx, index) 
		for( idx = 0 ; idx <= i ; idx++){
		 index = (idx + 1) * max_cols + (i + 1 - idx);
         input_itemsets[index]= maximum( input_itemsets[index-1-max_cols]+ referrence[index], 
			                             input_itemsets[index-1]         - penalty, 
									     input_itemsets[index-max_cols]  - penalty);

		}
	}
	printf("Processing bottom-right matrix\n");
    //Compute bottom-right matrix 
	for(i = max_cols - 4 ; i >= 0 ; i--){
       #pragma omp parallel for shared(input_itemsets, referrence, i,max_cols,penalty) private(idx, index) 
       for( idx = 0 ; idx <= i ; idx++){
	      index =  ( max_cols - idx - 2 ) * max_cols + idx + max_cols - i - 2 ;
		  input_itemsets[index]= maximum( input_itemsets[index-1-max_cols]+ referrence[index], 
			                              input_itemsets[index-1]         - penalty, 
									      input_itemsets[index-max_cols]  - penalty);
	      }

	}
    }
    CLOCK_End(chrono, omp_time);
   	

}
